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Abstract. We analyze spreading of a density front in the Kiintz-Lavallee model of 
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porous media. In contrast to previous studies, where unusual properties of the front 
were attributed to anomalous diffusion, we find that the front evolution is controlled 
by normal diffusion and hydrodynamic flow, the latter being responsible for apparent 
enhancement of the front propagation speed. Our finding suggests that results of 
several recent experiments on porous media, where anomalous diffusion was reported 
based on the density front propagation analysis, should be reconsidered to verify the 
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^ ! 1. Introduction 

The root-mean-square displacement, y(r 2 ), of a single diffusing particle at long times 
t is typically proportional to t a with a = 1/2. The same scaling law is also usually 
obeyed by the width of a region infiltrated by a diffusing front. However, there are 
also many natural phenomena where the exponent a differs from 1/2. These include, 
for example, water absorption in porous building materials [H [21 El H], copper sulfate 
diffusion into deionized water [51 E], transport in high polymer systems [7j, diffusion 
through synthetic membranes [8], and some aspects of surface diffusion [9]. Diffusion 
processes with a = 1/2 are ubiquitous in Nature and hence called 'normal', whereas 
those with a ^ 1/2 are quite rare and called 'anomalous'. 

There are several mechanisms that are known to bring about anomalous tracer [10] 
diffusion; for example, a power-law probability distribution of waiting times between 
individual steps of a random walk or a power-law probability distribution of distances 
covered by the diffusing particle at each jump [11] H2] . However, far less is known 
about anomalous chemical [10] diffusion, i.e. diffusion of macroscopic quantities of 
matter. Anomalous spreading of diffusion fronts is a many-body effect that results from 
complicated interactions between diffusing particles as well as between the particles and 
the surrounding (e.g. porous or fractal-like) medium. In many cases the nature of these 
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interactions remains unknown, so actually there is no satisfactory explanation of the 
anomalous diffusion in these systems [9]. Thus, any progress in the theory of anomalous 
chemical diffusion would be much welcomed. 

Recently Kiintz and Lavallee [H [131 El suggested that anomalous chemical diffusion 
can be a common phenomenon in systems with concentration-dependent diffusion 
coefficient of the diffusing entity. They based it on extensive computer simulations 
of a lattice-gas automata (LGA) model p21E|, and a refined analysis of a CUSO4 front 
diffusion into deionized water experiment [6]. Were this conjecture correct, it could be 
used to identify the cause of anomalous spreading of concentration fronts in many real 
systems. 

In this paper we focus on the Kiintz and Lavallee (KL) model and propose a much 
simpler interpretation of the simulation results obtained in [13j [14]. Our approach 
indicates that the diffusion in the KL model is normal. The key to the proper 
interpretation of the simulation results is the fact that the KL model is an extension 
of the Frish-Hasslacher-Pomeau (FHP) lattice gas model of fluid flow. Even though 
Kiintz and Lavallee applied it in a porous medium, which significantly reduced the 
flow velocity, hydrodynamic effects were not eliminated altogether. We found that the 
fluid flow velocity is large enough to enhance the effective speed of the diffusing front, 
apparently enlarging the value of the scaling exponent a above 1/2. We show that when 
the flow is taken into account, the exponent a assumes its expected, "normal" value 
1/2. 

In a broader perspective, our study shows that application of "hydrodynamic" 
lattice-gas automata models to study diffusion in porous media requires a great caution. 
Special attention must be paid to modelling the porosity of the medium to ensure that 
there is no macroscopic flow through the system. If this is impossible, the flow, even of 
minute magnitude, must be explicitly taken into account in the analysis of results. 

The structure of the paper is as follows. In Section [2] we briefly introduce the Kiintz- 
Lavallee model and the main results obtained thus far. In Section [3] we present our 
interpretation of these results, and show that actually there is no anomalous diffusion in 
the system. Section H] presents discussion of the major properties of the model. Finally, 
Section [5] is devoted to conclusions. 

2. The Kiintz-Lavallee model 

The Kiintz-Lavallee model [HI [131 El IS] describes diffusion of a fluid in a porous 
medium and is an extension of a deterministic Frish-Hasslacher-Pomeau (FHP 5 ) lattice- 
gas automata model of fluid flow [16]. The system is reduced to a two-dimensional 
triangular lattice, L x lattice units (l.u.) long and L y \/3/2 l.u. wide. The total number 
of lattice nodes is thus N = L x L y . They are occupied by indistinguishable particles, 
whose velocities belong to a discrete, 7-element set (Fig. CO). Any particles occupying 
the same lattice site must have different velocities, so that each node can be occupied 
by at most 7 particles. The time is discrete and at each time step particles either stay 
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at rest (if their speed is 0) or jump to the adjacent site pointed at by their velocity 
vector. Particles arriving simultaneously at the same node may collide. The collision 
rules for the FHP5 model are depicted in Fig. [H they conserve the mass and momentum. 
Additionally, a porous medium is modelled by assuming that a fraction c s of randomly 
chosen lattice nodes is permanently occupied by the so called specular scatterers (such 
scatterers behave like tiny rigid squares parallel to the "x" axis [15]). The system is 
confined between rigid, impenetrable walls along the "x" direction, whereas periodic 
boundary conditions are assumed along the "y" direction. A characteristic feature of 
the model thus constructed is a strong dependence of the diffusion coefficient D on the 
reduced local particle concentration c p3], which is defined as the average number of 
particles per node divided by 7. It turns out that D assumes a minimum value at c ~ 0.2 
and diverges to infinity as c goes either to or 1 [13J. 




Figure 1. Collision rules in the FHP5 model. Filled circles represent particles at rest. 

Initially the system is uniformly filled with particles so that the reduced 
concentration is Ci throughout the system. To this end 7c2-/V particles are distributed 
randomly in the system in such a way that their mean velocity vanishes and any particles 
occupying the same node have different velocities. Next, a narrow strip of nodes at the 
left-hand side of the system is uniformly filled with additional particles so as to increase 
the reduced concentration to c\ > C2. In this way a step-like nonuniform initial condition 
is formed, 

c(x, 0) = c\H{—x), x < 0, c(x, 0) = c 2 H(x), x > 0, (1) 

where H(x) is the Heviside step function. The origin of the reference system is set to 
the interface between the areas of high and low initial concentration, so that x = 
corresponds to the initial location of the step. 

During the simulation, by injecting new particles if necessary, the reduced 
concentration in the boundary strip x < is maintained fixed at c\. This corresponds 
to boundary conditions that are often used in studies on water infiltration into a porous 
medium: 

c(x, t) = c\ for x = 0, c(x, t) — > c 2 as x — > 00. (2) 
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Figure 2. Reduced concentration profile c(x, t) as a function of position x for 
t = 8000 x 2 fe time steps, k = 0,1,..., 6. The boundary conditions are Ci = 0.9 
and C2 = 0.2. 



3. Results 

We performed extensive simulations of a concentration front propagation using a 
triangular lattice with L x = 8000, L y = 200, the density of scatterers c s = 0.08, and 
the initial and boundary conditions determined by c\ = 0.9 and C2 = 0.2. The results 
for t = 8000 x 2 k , k = 0, 1, . . . , 6, are depicted in Figure [2] as a semi-logarithmic plot. 
Both the shape of the consecutive curves and the distances between them appear to be 
almost the same, which suggests that the concentration profile c(x, t) can actually be 
expressed function of a single variable x/t a . 

Validity of the scaling hypothesis can be verified by plotting several concentration 
profiles, obtained at different times, as functions of x/t a . If the scaling holds, all these 
plots should collapse into a single curve. Such analysis was already performed in [13J, 
where it was found that a ~ 0.55. This was interpreted as an evidence that the transport 
in the LGA model is super diffusive (a > 1/2). 

This finding, however, leads to a paradox: how addition of scatterers can enhance 
transport of particles? According to Kiintz and Lavallee [13] , this paradox can be 
explained by two hypotheses. The first one attributes superdiffusion of particles to a 
very strong dependence of their diffusivity on concentration. The second one says that 
anomalous diffusion is actually a transient phenomenon: in the limit of time t — > oo the 
scaling exponent a will extremely slowly decrease to 1/2. 

However, this argumentation neglects a well known fact, widely used in the 
Boltzmann-Matano method jTOj, [T7j, that any solution to the diffusion equation 

dc(x,t)_ d ( D{c) d±A\^ (3) 



dt dx \ dx 

with the initial and boundary conditions ([1]) and ([2]), satisfy the x/t 1 ^ 2 scaling for all 
t, irrespective of the form of D(c) [17]. Consequently, absence of this scaling indicates 
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Figure 3. The "x" component of the velocity, v x , as a function of the distance from 
the origin, x, for time t = 8000 x A k , k = 0, 1, 2. 



that equation ([3]) cannot give a proper description of the concentration front dynamics in 
the model. Besides diffusion, there must be another transport mechanism that controls 
particle infiltration into the low-concentration area. We conjecture that this is the fluid 
flow induced by a pressure gradient between the high- and low-concentration areas of 
the system. 

To check this hypothesis, we started from detailed analysis of spacial and temporal 
properties of the "x" component of the velocity vector, v x (x,t). Figure [3] depicts the 
velocity profile at several times t. It shows that the region x > can be split into 
two parts characterized by different behavior of v x . The first region coincides with the 
particle infiltration area. The system already "feels" the pressure gradient there and 
responds to it by developing a small but definitely nonvanishing flow towards the low 
concentration area. The second region is located further away from the infiltration area. 
The system remains there close to the initial state, with v x fluctuating around 0. 

In order to study the concentration front dynamics, one thus have to take into 
account the front velocity, V{. In general, this quantity is both x- and t-dependent. 
However, we are interested only in its approximate value at time t, and estimate it as a 
mean weighted with c(x,t) — 

POO 

/ v x (x,t)[c(x,t) - c 2 ] dx 

Vt® = ^—755 • (4) 

/ [c(x,t) - c 2 ] dx 
Jo 

In this definition we utilize the fact that the weight function quickly vanishes outside 
the front region (Fig. [2]), while the averaged quantity, v x , remains almost constant there 
(Fig. Ej). Moreover, the integral form of the definition helps to significantly reduce the 
influence of relatively large statistical noise present in the data for v x . 
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Figure 4. The mean front velocity, V{, as a function of time. The dashed line is a 
guide to the eye with a slope —1/2. 



Figure H] shows that the initial magnitude of Vf is very high, almost reaching the 
maximum possible value 1. At larger times V{ decreases as with /3 < 1/2 very weakly 
depending on t. Consequently, the distance the front moves due to the flow, 

s{t) = [ v f (r) dr, (5) 



for large t grows as i 1- ^, i.e. a bit faster than the typical diffusive length m t 1 / 2 . This 
suggests that the front dynamics is controlled by both diffusion and flow. Actually, 
our simulations showed that after t = 256000 time steps the front displacement was 
Ax « 3750 l.u. (see Figure [2]), of which s(t) ~ 1440 l.u, that is, about 38%, was caused 
by the fluid flow. 

To separate diffusion from flow, we investigated the front propagation using a mobile 
reference frame (x',t), with 

x' = x-s(t). (6) 

This new coordinate system moves along the "x" axis with velocity Vf(t), so as to 
minimize effects of the fluid flow. Assuming that in this new reference system the 
front dynamics is dominated by (normal) diffusion, the concentration profiles should 
asymptotically obey a similarity relation 

c(x' : t) = f((x'-x )/Vt), (7) 

where / is a similarity function and xq is a constant. To verify this conjecture, we first 
used the data for c(x', t) = 0.5 to estimate Xq ~ —190, and then plotted c as a function of 
(x' — Xo)/\rt for several times t. The results, obtained for 8000 <t< 512000, are shown 
in Figure [5j They confirm that the concentration profiles asymptotically satisfy ([7]), as 
for t > 64000 the scaled-up profiles are practically indistinguishable. 
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Hydrodynamical flow can also explain another feature of the model — discontinuity 
of fluid concentration at the boundary x — 0. Even though the boundary condition (J2J) 
forces the concentration to be fixed at c\ for x < 0, its values at x > are significantly 
smaller than c\ (see Figure (j2]), obtained for c\ = 0.9). Discontinuity in c, Ac, is 
related to discontinuity in the flow velocity Av ~ Vi (v x = for x < and v x ~ Vi 
for x > 0). Let c + = ci — Ac denote the concentration at x = 1 l.u. The rate of 
diffusive particle transfer per unit length through the interface x = is approximately 
equal to 3ci — 3c + = 3Ac (the factor 3 is the number of velocity directions in the LGA 
model that transfer particles through the interface). This must be balanced by particle 
current caused by the fluid flow, which can be approximated by 7c + V{ (the factor 7 is 
the number of possible particle velocities). Thus, 3ci — 3c + « 7c + V{, or 

Ac ^ 8 
3 + 7vf 

Comparison of this formula with the actual data obtained in simulations is shown in 
Figure (Q. The agreement is fairly good. 



4. Discussion 



Several features of the model deserve closer attention. First, the model employs specular 
scatterers to mimic a porous medium as well as to reduce hydrodynamic effects [13J. 
However, the effect of such scatterers on the diffusion along the "x" axis is very peculiar: 
the "x" component of the momentum is changed only for those particles, whose velocity 
vector is parallel to the "x" axis. For this reason specular scatterers are widely used 
in LGA models to implement slip-free walls parallel to the x axis. However, in the KL 
model the scatterers are distributed randomly with a relatively small concentration to 
form isolated, point-like islands. Such scatterers are inefficient in breaking temporal 
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Figure 6. Reduced concentration drop, Ac, at the boundary x = as a function of 
the front velocity Vf, obtained from simulations (□) and equation ([5]) (solid line). 



correlations in the "x" component of momentum of particles that hit them. Thus, on 
the one hand, randomly distributed specular scatterers have a very limited impact on 
the diffusivity along the "x" axis. On the other hand, the mean-free path, and hence 
the diffusion coefficient in the original, scatterer-free FHP model diverge to infinity as 
the reduced particle concentration goes to or 1 [13J. Combination of these two effects 
explains the fundamental property of the KL model: extremely strong dependence of 
the diffusion coefficient on concentration [T3] . 

Second, the present model is based on the FHP model, where particle-particle 
collisions preserve mass and momentum to mimic fluid flow. Local conservation of 
momentum leads to hydrodynamic flow from high-pressure (high-concentration) to low- 
pressure (low-concentration) areas. Therefore, hydrodynamical flow in any FHP-based 
model of porous media is inevitable and its effects should be always carefully taken into 
account, especially when specular scatterers are used to model a porous medium. 

Third, assume that as soon as the flow velocity in the front becomes sufficiently 
small, the particle flux q associated with it obeys Darcy's law [T8] 



where k is the permeability of the medium, jj, denotes the fluid viscosity, and VP is 
the pressure gradient. Taking into account that q oc v x [18J, /i oc D^ 1 (Stokes-Einstein 
relation), and P oc c [15] . one arrives at 



If D is very large, even a very small concentration gradient can cause a significant 
fluid flow. Therefore, one can expect strong hydrodynamical effects to appear in any 
hydrodynamical model of porous medium with diffusion coefficient D(c) diverging to 
infinity at some value(s) of c. 




(9) 



v x oc DVc. 



(10) 
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Fourth, the value of |x | in the similarity relation ([7]) turned out quite large. 
Similarly, the relaxation time to the asymptotic, long-time regime in which the similarity 
holds, r re i ax , is also large. These two effects may be related to the fact that in all FHP 
models v x , which constitutes the left-hand side of relation (flOj) . is limited from above. 
However, the right-hand side of (TTUl) is a product of two quantities that at early stages of 
the front propagation assume very large values in the front region. Thus, the system can 
obey Darcy law only after a long relaxation time necessary to reduce the concentration 
gradient Vc to a value of order 1/D. This explains why in the previous studies the 
front dynamics was found "anomalous" only if both c\ was close to 1 and the specular 
scatterers were used P, [131 Ej- If an y of these conditions was not satisfied, the value 
of diffusion coefficient D took on relatively small values in the front region and the 
relaxation time was small enough to allow the simulations to reach the asymptotic 
regime. Note also that large values of \xq\ and r re i ax are in accord with several recent 
experiments on sorptivity of building materials, where a short-time deviation from the 
t 1 ! 2 behavior was found (see [I] and references therein). 

Fifth, the actual front velocity at early times turned out smaller than its value 
extrapolated from the long-time, "Darcy" asymptotics (see Fig. HI). As a consequence, 
Xq is negative and the front dynamics at early times appears "superdiffusive" (a > 1/2). 
In systems where the early-time velocity is larger than the short-time extrapolation of 
the Darcy asymptotics, xq would be positive and the early-time front dynamics would 
appear "sub diffusive" (a < 1/2). 

Finally, Figure H] suggests that at large times the mean front velocity, Vf, becomes 
proportional to l/y/i, which yields a characteristic length scale J^ =0 l/-y/rdT oc \ft. 
This implies that the characteristic lengths associated to both diffusion and fluid flow 
asymptotically scale with time in the same way, as \fi. Consequently, the KL model 
is an example of a system where measurements based only on analysis of asymptotic 
concentration profiles cannot differentiate between diffusive and hydrodynamic effects. 

5. Conclusions 

We have shown that the unusual dynamics of the concentration front in the KL model 
of a porous medium can be fully explained as a combined effect of hydrodynamic flow 
and normal diffusion. Our argumentation is not only simpler than that proposed by 
Kiintz and Lavallee [H [131 E] > wno used the concept of anomalous diffusion, but also 
allows to explain more physical properties of the model, e.g. a concentration drop at the 
boundary x — 0. 

From our analysis the following picture of the concentration front dynamics in 
the model emerges. A concentration gradient across the front leads to a pressure 
gradient which, in turn, forces the FHP fluid to flow. This flow is attenuated by 
scatterers that mimic a porous medium, so that hydrodynamical effects in most cases are 
negligible, especially at large times. However, randomly distributed specular scatterers 
used in the KL model are very inefficient at reduced concentrations close to 1. In this 
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particular case the viscosity goes to 0, and both diffusivity and mobility along the "x" 
axis tend to infinity. Under these conditions, hydrodynamical flow significantly affects 
the front dynamics at all times. In particular, it renders the relaxation time to the 
asymptotic regime, where Darcy law holds, very large, and significantly changes the 
position of the Boltzmann-Matano interface Xq. The front dynamics is controlled by 
both diffusion and fluid flow, and at very large times the characteristic length scales 
of both processes become proportional to \/t. Consequently, at very large times the 
concentration profile can be expressed as a function of a single variable x/y/i, which 
could easily be misinterpreted that the front dynamics is governed only by diffusion. 

While the KL model cannot be regarded as a model of anomalous diffusion, it 
remains an interesting model of fluid flow in a porous medium. Its main advantage is 
ability to reproduce some non-trivial properties of concentration profiles found in several 
recent experiments. Similarity between the model an experimental results suggests that 
the anomalous front behavior observed e.g. in experiments on building materials, may 
have nothing to do with anomalous diffusion, but is caused by some hydrodynamic 
effects that accompany normal diffusion. 
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